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Abstract 



Consider a linear regression model with regression parameter /3 = (Pi, . . . , /3 P ) and 
independent normal errors. Suppose the parameter of interest is 9 = a T (3, where a 
is specified. Define the s- dimensional parameter vector r = C T f3 — t, where C and t 
are specified. Suppose that we carry out a preliminary F test of the null hypothesis 
H : r = against the alternative hypothesis Hi : r ^ 0. It is common statistical 
practice to then construct a confidence interval for 9 with nominal coverage 1 — a, 
using the same data, based on the assumption that the selected model had been given 
to us a priori (as the true model). We call this the naive 1 — a confidence interval 
for 9. This assumption is false and it may lead to this confidence interval having 
minimum coverage probability far below 1 — a, making it completely inadequate. 
Our aim is to compute this minimum coverage probability. It is straightforward 
to find an expression for the coverage probability of this confidence interval that 
is a multiple integral of dimension s + 1. However, we derive a new elegant and 
computationally-convenient formula for this coverage probability. For s = 2 this 
formula is a sum of a triple and a double integral and for all s > 2 this formula 
is a sum of a quadruple and a double integral. This makes it easy to compute 
the minimum coverage probability of the naive confidence interval, irrespective of 
how large s is. A very important practical application of this formula is to the 
analysis of covariance. In this context, r can be defined so that Hq expresses the 
hypothesis of "parallelism". Applied statisticians commonly recommend carrying 
out a preliminary F test of this hypothesis. We illustrate the application of our 
formula with a real-life analysis of covariance data set and a preliminary F test 
for "parallelism". We show that the naive 0.95 confidence interval has minimum 
coverage probability 0.0846, showing that it is completely inadequate. 

Keywords: analysis of covariance, naive confidence interval, preliminary F test, test 
for parallelism. 
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1. Introduction 

Consider the linear regression model Y = X(3 + e, where Y is a random n- vector of 
responses, X is a known n x p matrix with linearly independent columns, (3 is an 
unknown parameter p-vector and e ~ N(0, <7 2 I n ) where a 2 is an unknown positive 
parameter. Suppose that the parameter of interest is 9 = a T (3 where a is a given 
p- vector (a ^ 0). We seek a 1 — a confidence interval for 6. 

Let the s-dimensional parameter vector r be defined to be C T (3 — t where C 
is a specified p x s matrix (s < p) with linearly independent columns and t is a 
specified s-vector. Suppose that a does not belong to the linear subspace spanned 
by the columns of C. Also suppose that we carry out a preliminary F test of the 
null hypothesis H : r = against the alternative hypothesis Hi : r ^ 0. It is then 
common statistical practice to construct a confidence interval for 6 with nominal 
coverage 1 — a, using the same data, based on the assumption that the selected 
model had been given to us a priori (as the true model). We call this the naive 
1 — a confidence interval for 6. In Section 2, we provide a convenient description of 
this confidence interval. This assumption is false and it can lead to the naive 1 — a 
confidence interval having minimum coverage probability far below 1 — a, making it 
completely inadequate. Our aim is to compute this minimum coverage probability. 
For s — 1, the preliminary F test is equivalent to a t test. The case of a single 
preliminary t test has been dealt with by Kabaila and Giri (2009, Theorem 3). So, 
in the present paper, we restrict attention to the case that s > 1. 

Straightforward application of the methodology of Farchione (PhD thesis, 2009, 
Section 5.7) leads to an expression for the coverage probability of the naive 1 — a 
confidence interval, for a given value of an s-dimensional parameter vector, that is 
a multiple integral of dimension s + 1. Finding the minimum coverage probability 
using this formula becomes increasingly cumbersome as s increases due to both the 
need to (a) evaluate multiple integrals of dimension s + 1 and (b) the need to search 
for the minimum over a space of dimension s. 

In Section 3, by a careful consideration of the geometry of the situation, we derive 
a new elegant and computationally-convenient formula for the coverage probability 
of this confidence interval for given parameter values. For s = 2 this formula is 
a sum of a triple and a double integral and for all s > 2 this formula is a sum 
of a quadruple and a double integral. This formula also shows that the coverage 
probability is a function of a two-dimensional parameter vector, irrespective of how 
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large s is. This makes it easy to compute the minimum coverage probability of the 
naive confidence interval, irrespective of how large s is. Another important aspect 
of this formula is that it can be used to delineate general categories of a, C and X 
for which the naive confidence interval has poor coverage properties. 

A very important practical application of this formula is to the analysis of co- 
variance. In this context, r can be defined so that H expresses the null hypothesis 
of "parallelism". In the applied statistics literature on the analysis of covariance 
it is commonly recommended that a preliminary F test of the null hypothesis of 
"parallelism" be carried out. See, for example, Kuehl (2000, p. 563), Milliken and 
Johnson (2002, pp. 14 - 17) and Freund et al (2006, pp. 363 - 368). For an analysis 
of covariance, we can choose a so that the parameter 9 is the difference in expected 
responses for two specified treatments, for the same specified values of the covariates. 

In Section 4, we illustrate the application of the results of the paper with a real- 
life analysis of covariance data set and a preliminary F test for "parallelism" . We de- 
fine #tobe (expected response to treatment 1) — (expected response to treatment 2), 
evaluated at the same specified value of the covariate. We show that the naive 0.95 
confidence interval for 9 has minimum coverage probability 0.0846, for this speci- 
fied value of the covariate. This shows that this confidence interval is completely 
inadequate, for this specified value of the covariate. 

2. Description of the naive confidence interval 

In this section we provide a convenient description of the naive 1 — a confidence 
interval constructed after the preliminary F test. Let j3 denote the least squares 
estimator of f3. Define R(f3) = (Y — Xf3) T (Y - X/3). Let m = n - p. Define 
£2 = R0)/ m = (Y - Xf3) T {Y - X$)/m. Also, define Q = a T (3 and f = C T /3 - t. 
We suppose that the columns of the matrix C are linearly independent. We also 
suppose that a does not belong to the linear subspace spanned by the columns of 
C. Now define the (s + 1) x (s + 1) matrix 



Note that v u = a T (X T X)- 1 a, v 2 i = C T (X T X)- 1 a and V 22 = C T {X T X)- 1 C. 



Define (3* to be the value of j3 minimizing R(j3) subject to the restriction that 



V 
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t = C T f3 -t = 0. As is well known (see e.g. Graybill, 1976, p.222) 



{3* = - (X T X)- 1 C(C T (X T X)- 1 C)- 1 r 
R0*) = R0) + 0- (3*) T (X T X)0 - n- 

The standard test statistic for testing H : r = against Hi : r ^ is 

p 0- P*) T X T X0 - p*)/s r T V^f/s 
t 2 t 2 

This test statistic has an F StTn distribution under H . Suppose that we reject H 

when F > £ and accept H otherwise, where £ is a specified positive number. 

Define 6* = a T /3*. Also define the quantile t(m) by the requirement that P( — 

t(m) < T < t(m)) = 1 - a for T ~ t m . The naive 1 — a confidence interval for 8 is 

obtained as follows. 

Suppose that F > I. The confidence interval is constructed on the assumption 
that r = is not necessarily true. In this case, the naive 1 — a confidence interval 
is the usual 1 — a confidence interval for 9 based on fitting the full model, 

/= [Q - t{m)y/v^% Q + t(m)y/v^t]. 

Now suppose that F < £. The confidence interval is constructed on the assump- 
tion that t = 0. If t = then R{f3*)/a 2 ~ xl +s and Var(0*) = a 2 (v ll -v% l V 22 1 v 2 i) ■ 
Note that 6* and R(/3*) are independent random variables. We use the notation 
[a ± b] for the interval [a — b, a + b] (b > 0). In this case, the naive 1 — a confidence 
interval for 9 is 



J 



Q*±t(m + s) 
@*±t(m + s) 



R(/3* 



m 



vn - vZ 1 V 22 1 v 2 i 



— — \hu - v{ x V 22 v 21 



3. The coverage probability of the naive confidence interval 



Define b 



-l/2 T ,-l/2 



V 22 v 2 \ and W = 11/ a. Let fw denote the probability density 



function of W. Define 



+ b 2 . Thus 



a T (X T X)- 1 C(C T (X T X)- 1 Cy 1 C T (X T X)- 1 a 
a T {X T X)- 1 a ' 
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Since Var(0*) = cr 2 [vu — f^V^ 1 ^!.) > 0, ||6|| G [0,1]. The assumption that the 
vector a does not belong to the linear subspace spanned by the columns of C implies 
that H&ll > 0. So, we may assume that ||6|| G (0, 1]. Now define 

i(x,w; \\b\\) — P( — t(m)w + x < Z < t(m)w + x), 

where Z ~ N(0, 1 - ||6|| 2 ), and 

j (x, y, «,; ||6||) = P ( x - t(m + s) J™* ' + V y/l - \ \b\ \ 2 < Z < 
\ V m + s 



Imw 2 + y r — 
x + t(m + s)\ a/ 1 - 

' " m + s 



Define /r to be the probability density function of vP? when P 2 ~ x 2 . Let B(a, b) 
denote the beta function. Define the probability density function fa to be 

7T 



^-sin- 2 ( 7 rt 1 ) far0<t!<l 

/Ti(ClJ = < ^2' 2 J 

I otherwise. 
For s > 3, define the probability density function fa 2 to be 

71 -,«-3, 



, ^ j mi s^ sin (7ri2) for0 ^^<i 

I otherwise. 

Let 7 = (l/o")V r 22 1 ^ 2 r. Define /q to be the probability density function of a noncen- 
tral chi squared distribution with s degrees of freedom and noncentrality parameter 
I 2 . Also define 



d(h,r;s, ||7| 



|7|| 2 + 2||7||rcos(27rt 1 ) + r 2 for s = 2 
|7|| 2 + 2||7||rcos(7rfi) + r 2 for s > 3. 



Define the unit vector u b = (1/||6||)6. When ||7|| > 0, define w 7 = (I/H7ID7 and 
then ip = ujiiry. Also define ip = 1 when ||7|| = 0. Now, when ||7|| > 0, define 

kfaip) = ^cos(27rti) + y/l - i) 2 sin(27rti) 

^ ^_ ^ I if) cos(7rti) + y/l — ip 2 sin(7rti) cos(27rt 2 ) for s = 3 
I ip cos(7rti) + y/l — ip 2 sin(7rti) cos(7rt 2 ) for s > 4. 

The following is the main result of the paper. 
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Theorem 1. TTie coverage probability of the naive 1 — a confidence interval for 6 
is P{6 G I,F > £) + P{6 G J, F < £) . A computationally- convenient expression for 
the second term in this sum is 

P(6eJ,F<£)= / / j(\\b\\MiP,q,w;\\b\\)f Q (q)f w (w)dqdw (2) 
Jo Jo 

and computationally- convenient expressions for P{6 G I,F > £) are as follows. Let 
u = ^d(ti,r; s, ||7||)/s£. For s = 2, 

pi poo pu 

P(9eI,F>£)= / / i(-\\b\\rk(t 1 -i;),w;\\b\\)f w (w)f R (r)dwdrdt 1 (3) 
Jo Jo Jo 

Fors>3 and \ \-f\\ > 0, P(9 G I,F > I) is equal to 

^1 pi poo pu 

i(-\\b\\r k(t u t 2 ; s, ip),w; \ \b\\) f w {w) f R (r) f Tl (ti) f T2 (t 2 ) dw dr dt x dt 2 

(4) 

For s > 3, 1 17| | > and ip G {-1, 1} ; 

pi poo pu 

P(9eI,F>£)= / / / i(-\\b\\rcos(7Tt 1 ),w;\\b\\)f w (w)f R (r)f Tl (t 1 )dwdrdt 1 
Jo Jo Jo 

For s > 3 and \ \ = 0, 

c-l poo pr/\/~sl 

P(8eI,F>£)= / / l (-\\b\\rco S (2nt 1 ),w;\\b\\)f w (w)f R (r)dwdrdt 1 
Jo Jo Jo 

Note that for given \ \b\\ (which is determined by a, C and X ) and m, s, £ and a, the 
coverage probability of the naive 1 — a confidence interval is a function of (| |7||, '0) ■ 

The proof of this theorem is presented in Appendix A. 

The formulas given in Theorem 1 have three attractive features. The first of 
these is that, irrespective of how large s is, these formulas involve, at most, a 4- 
dimensional integral. The second of these features is that the numerical evaluation 
of these integrals, reviewed in Appendix B, is very straightforward. Thirdly, for 
given m, s, £, a and ||6||, the coverage probability of the naive confidence interval is 
a function of the two-dimensional parameter vector (||7||,V ; )- These features make 
it is easy to compute the minimum coverage probability of the naive 1 — a confidence 
interval for given to, s, £, a and ||6||. Finally, Theorem 1 has the following corollary 

Corollary 1. For given m, s, £ and a, the infimum over (IHI,^) of the coverage 
probability of the naive 1 — a confidence interval for 6 is a function of 
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For the numerical example described in the next section, m = 4, s = 4, £ = 
6.5914 (corresponding to a 0.05 level of significance of the preliminary F test) and 
a = 0.05. For these values of m, s and £, the minimum coverage probability of 
the naive 0.95 confidence interval, as a function of ||6|| G (0,1], is as shown in 
Figure 1. All of the computations presented in the present paper were performed 
with programs written in MATLAB, using the optimization and statistics toolboxes. 
We note the decrease in the minimum coverage probability of this naive confidence 
interval with increasing ||6||. We can use Corollary 1 to delineate general categories 
of a, C and X (via their relationship to ||6||) for which this naive 0.95 confidence 
interval has poor coverage properties. Specifically, this naive confidence interval will 
have poor coverage properties for those values of a, C and X such that 

a T {X T X)- 1 C(C T {X T X)- 1 C)~ 1 C T {X T X)- 1 a 
a T (X T X)- l a 

is greater than 0.7, say. 





0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Hb|| 

Figure 1: Plot of the minimum coverage probability of the naive 0.95 confidence 
interval, as a function of ||6|| e (0, 1], for m = 4, s = 4, £ = 6.5914 and a = 0.05. 
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4. Application to a real-life data set 



In this section we consider the real-life analysis of covariance data set due to Chin et 
al (1994) and analysed by Yandell (1997, Chapter 17), who makes this data available 



at the website http://www.stat.wisc.edu/~yandell/pda/. This data is listed in the 
Table 1. It consists of the observed response (weight gain) for a given treatment 
and value of the covariate (feed intake). There are 4 possible treatments, numbered 
1 to 4. 



treatment 


weight gain 


feed intake 


1 


1416.1 


2451.75 


1 


1447.0 


2546.00 


1 


1509.6 


2657.00 


2 


1497.8 


2452.10 


2 


1469.9 


2404.90 


2 


1469.4 


2479.90 


3 


1510.1 


2788.70 


3 


1423.0 


2655.50 


3 


1295.9 


2366.40 


4 


1354.8 


2578.80 


4 


1326.8 


2384.40 


4 


1335.1 


2477.50 



Table 1 : The observed response (weight gain) for a given treatment and value of the 
covariate (feed intake). Source: j http: / /www.stat.wisc.edu/ ~yandell/pda/[ 

We use the following linear regression model for this data: 

Yij /ij -|- fii{xij X..) -\- Ey 

where is the response of the j th experimental unit (j = 1, . . . , 3) that is receiving 
the i th treatment (i = 1, ... ,4), when the covariate takes the value x^. The 
are independent and identically N(0, a 2 ) distributed and a 1 is an unknown positive 



parameter. The fii and /3, (i 



,4) are unknown parameters. Also, x.. denotes 



the mean of the Xij (i 



1,...,4; j = l,...,3). 
We express this model in the form Y = X/3 + e, where 

Y — (Yn, Y u , Yis, Y 2 i, Y 2 2, Y 2 3, Y 3 i, Y 32 , Y 33 , y 41 , y 42 , 143) 
£ — (^ll) £12, £i3) £ 2i, £22, £23, £31, £32, £33, £41, £42, £43) , 



j3 = (/ii, /i 2 , /X3, /U 4 , fix, $2, $3, M and X is the obvious 12 x 8 matrix. 
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As considered by Yandell (1997, p. 271), we carry out a preliminary test of the 
null hypothesis H : f3\ = (3 2 = /?3 = $a against the alternative hypothesis Ha that 
the $i are not all the same, using an F test. Suppose that we use a 0.05 level of 
significance for this preliminary test. We express H as r = and H A as r ^ 0, 
where r = C T (5 and 



Define the parameter of interest 9 as follows. Let Y£ and Y 2 * denote the responses 
of two experimental units, receiving treatments 1 and 2 respectively, for the same 
value x* of the covariate. In other words, 



where e\ and e 2 are independent and identically N(0, a 2 ) distributed random vari- 
ables. Let 6 = E(Y?) - E(Y 2 *) = ^ - [i 2 + 0i - p 2 )(x* - x..). Thus, 9 = a T (5 where 
a = (l, -1,0,0, (x*-x..), -(x* -x..),0, 0). We suppose that x* — x.. = 125.39, which 
is the maximum value of \xa — x..\ for the data. 



For this situation, ||6|| = 0.96869 and so the minimum coverage probability of 

the naive 0.95 confidence interval is 0.0846. This shows that this confidence interval 
is completely inadequate, in this situation. 

5. Discussion 

The poor coverage properties of naive confidence intervals found in this paper are 
presaged by the poor coverage properties of naive confidence intervals found in the 
context of a preliminary best subset variable selection by minimizing an AlC-type 
criterion, see e.g. Kabaila (2005), Kabaila & Leeb (2006) and Kabaila & Giri (2009). 
Apart from the form of preliminary model selection used, minimum AIC versus an F 
test, these papers differ from the present paper in that the present paper provides a 
method for computing the minimum coverage probability, whereas Kabaila (2005), 
Kabaila & Leeb (2006) and Kabaila & Giri (2009) provide only upper bounds on 
the minimum coverage probability of the naive confidence interval. 



- 1 1 
C T =0000 -1010 
- 1 1 



Y* = m + f3 l (x 

y; = ii 2 + $ 2 {x 



x..)+e\ 
x..)+e* 2 
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Appendix A: Proof of Theorem 1 

In this appendix, we prove Theorem 1. Define G = (0 — 9) / {a-sjvu) and H 
(l/cr)V 2 2 t. Let fn denote the probability density function of H. Note that 



G 
H 



N 



1 b T 
b L 



(5) 



where I s denotes the s x s identity matrix. Thus the distribution of G, conditional 
on H = h, is N(b T (h - 7), 1 - ||6|| 2 ). Note that (G,H T ) T and W are independent 
random vectors. We use the notation 



1 if A is true 
if A is false 



where A is an arbitrary statement. This is similar to the Iverson bracket notation 
(Knuth, 1992). 

By the law of total probability, the coverage probability of the naive 1 — a 
confidence interval is 

P(0 eI,F>£) + P(6eJ,F<£). 

We divide the remainder of the proof into 2 parts. 
Part 1: expression for P(0 £ I,F > £) 

Suppose that ||7|| > 0. We prove the validity of the expressions ()3]) and @ for 
P(0 6 I,F > £). The proofs of the validity of the other expressions for P{9 G 
I,F > C.) (given in Theorem 1) are similar and are omitted, for the sake of brevity. 
Now {9 e 1} = {-t(m)W <G< t(m)W} and F = H T H/(sW 2 ). Thus 

P(9 eI,F>£) 

/ H T H \ 

= P -t(m)W <G< t(m)W, > £ 

\ sW 1 J 

(H T H \ 
-t{m)W <G< t(m)W, —— >£ W = w, H = h) f w (w) dw f H (h) dh 
sW 2 J 

= J ■■ J J p (-t(m)w <G <t{m)w H = X {h T h > s £w 2 ) f w (w) dw f H (h) dh 

(6) 

Note that 

P(-t(m)w <G< t{m)w \H = h) 

= P( - t(m)w - b T (h - 7) < Z < t(m)w - b T (h - 7)) where Z ~ N(0, 1 
= i( - b T (h->y),w; 
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Thus 

/p poo 
...J j i(-b T (h- 1 ),w;\\b\\) X (h T h> s£w 2 ) f w (w)dwf H (h)dh 

= e(i(- b T {H - 7), W; \\b\\) X (H T H > siW 2 )) 

= e(i(- b T (H - j),W; \\b\\) X (W < y/HTH/sl)). (7) 

We now find a simple formula for this expected value. Since H ~ iV(7, J s ), 
H = 7 + RU where R is a nonnegative random variable and U is a random s- vector 
with the following distributions. The random vectors R 2 and U are independent, 
with R 2 ~ \ 2 S and U is distributed uniformly on the surface of the unit sphere in W. 
Thus6 T (#-7) =R\\b\\L b , where L b = u\ r U. Also, = ||7|| 2 + 2||7|| J RL 7 + J R 2 , 

where L 7 = « 7 [/. Hence 

© = - \\b\\RL b , W- \\b\\) X (W < ^(\\ 1 \\ 2 + 2\\ 1 \\RL 1 + R 2 )/stj). 

Note that i? and are independent random vectors. Define the ran- 

dom vector (Ti,T 2 ) to be such that T\, T 2 , R and W are independent and 2\ and 
T 2 have the probability density functions Jt x and /t 2 respectively, defined in Sec- 
tion 3. Define the unit s-vectors e 7 and e b as follows. The vector e 7 has 1 as 
its first component and zeros for the remaining components. The vector e b has 



first component ip, second component \f\ — ifj 2 and zeros for the remaining compo- 
nents. Because U is distributed uniformly on the surface of the unit sphere in W, 
(L 7 , Lft) has the same distribution as (e^U,e b U). It follows from Fang and Wang 
(1994, p. 49, pp. 306-306 and p. 308) that (Ui,U 2 ) has the same distribution as (a) 
(cos(27rTi), sin(27rTi)) for s = 2, (b) (cos(7rTi), sin(7rTi) cos(27rT 2 )) for s = 3 and 
(cos(7rTi), sin(7rT 1 ) cos(7rT 2 )) for s > 3. Thus, (L 7 , L b ) has the same distribution as 

(coB(2itTi), ^cos(2ttT 1 ) + yjl - ip 2 sin(27rTi) 

for s = 2, 

^cos(7rTi), ip cos(7rTi) + a/1 — ip 2 sin(7rTi) cos(27rT 2 ) 

for s = 3 and 



cos(7rTi), ?/> cos^T 1 !) + a/1 — ^ 2 sin(7rT 1 ) cos(7rT 2 ) 
for s > 3. Thus P(0 G J, F > £) is 

£?(i(-||6|| J RA;(T 1 ;^) J W;||6||)x(^< ^(^,^2, || 7 ||)/s^) 
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for s = 2 and 

E(i{ - \ \b\ \ Rk(T 1 ,T 2 ] s, V), W; \\b\\) X {W < \/ (d(T\, R; s,\\l\\)/s( 

for s > 3. This leads to the expressions and (jlj) for P(6 G I,F > £) given in the 
theorem. 

Part 2: expression for P(6 £ J, F < £) 

The derivation of the expression for P(8 G J, F < £) is based on ([T]) and the fact 



that O*, f and R((3) are independent random vectors. Define Q = (l/o -2 )^^ 1 
and note that Q has a noncentral chi squared distribution with s degrees of freedom 
and noncentrality parameter ||7|| 2 . Note that 



T 



{6 e,J} = <Z e 



,i ,mW 2 + Q r- 

6 T 7 ± t(m + s)\ ~ , — Vl 



m + s 

where Z, Q and are independent random variables and Z ~ -/V(0, 1 — ||fr|| 2 )- Also 

r f^f/a Q/a 



W 2 



Thus 

P(6 eJ,F< 



P z e 



T , mW 2 + Q r- 
b j±t{m + s)\ a/ 1 



o jo 



m + s 

j{\\ b \ \ IItIIV^w; \ \b\\) f Q (q) f w (w) dqdw 



by a method similar to that used in Part 1. 

Appendix B: Numerical evaluation of the integrals in 
Theorem 1 

We evaluate the integrals (j2J) and ([3]) in the statement of Theorem 1 as follows. We 
approximate fl2j) by 

"Ci ps£w 2 



./0 



J (I HI INI^w; INI) /<?(?) fw(w) dqdw 



(9) 



for an appropriately chosen value of c\. We bound the error of this approximation 
as follows. Since 1 1-7| | -0, q, w; \ \b\\) is a probability, 

"CO rsiw' 2 



< 



< 



< 



fw(w)dqdw 

fw{w) dw 
P(M > m c 2 ) where M ~ x 



ci JO 

oo 



Cl 
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We choose C\ sufficiently large that the right hand side is less than, say, 10~ 5 . 

To evaluate ([9]), we transform the region of integration to a rectangle as follows. 
Change the variable of integration q in (J9]) to q* = qj (s£w 2 ), so that (jHJ) is equal to 

r j(||6|||| 7 ||V,s^V,^l|6||) f Q {stw 2 q*)stw 2 f w (w)dq*dw. 



'0 JO 

The integrand is a smooth function of (q*,w) G [0, 1] x [0,ci] and so this double 
integral is easily evaluated by numerical integration. 
We approximate (j3J by 
'1 rc2 ru 



ri pc2 ru 

I I I t( - 

Jo Jo Jo 



rkfaip)^; \ \b\\) f w (w) f R (r) dwdrdti (10) 



for an appropriately chosen value of c 2 . We bound the error of this approximation 
as follows. Since i( — \\b\\ r k(ti; ift), w; \\b\\) is a probability, 

-1 poo ru 



rl roc ru 

< © - (HDD < / / / f w (w)dwf R (r)drf T (t)dt 

JO Jc 2 JO 

roo 

< / fn{r)dr 

J C2 

= P(R 2 > c 2 2 ) where R 2 ~ xl 



We choose c 2 sufficiently large that the right hand side is less than, say, 10~ 5 . 

To evaluate (fTU]) . we transform the region of integration to a rectangle as follows. 
Change the variable of integration w in (JTDj) to w* — w/ ^/d(t 1: r; s, | pyJU/sZ, so that 
( TTU|) is equal to 

r \( - \\b\\ r k(h;i>), v^^HDT^*; INI) fw(Vd(t^MWtw*) 

^/d(ti,r;s, \\^\\)/s£ f R (r) dw* dr dti. 



o Jo Jo 



The integrand is a smooth function of (w*,r,ti) G [0, 1] x [0,c 2 ] x [0, 1] and so this 
triple integral is easily evaluated by numerical integration. 

The evaluation of the other integrals in the statement of Theorem 1 is sim- 
ilar to the evaluation of the integrals (J2J) and (J3J). The evaluation of (j4]) using 
MATLAB requires special comment. In MATLAB, the highest dimensional inte- 
gral that one can evaluate using a built-in MATLAB function is a triple integral. 
We evaluate the quadruple integral 01]) using MATLAB as follows. As before, let 
u = a/g^x, r; s, \\j\\)Jsi. Define 

g(w*,tx,r)= / i( - \\b\\ rk(t 1 ,t 2 ; s, ip), uw*\ \\b\\) f Ta (h) dt 2 
Jo 
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The integrand on the right-hand-side is a very smooth function of t^- We evaluate 
g(w*,ti,r), to a good approximation, using a compound Simpson's rule with a 
specified number of subdivisions of the interval [0, 1]. The quadruple integral (@J is 
approximated by 

pc 2 pi pi 

/ / / g{w\t u r) f w (uw*)uf Tl (ti) f R (r)dw* dt x dr 
Jo Jo Jo 

which is evaluated using the MATLAB built-in function triplequad. 
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